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C^ • We study the planar motion of test particles in gravitational fields produced 

^^ ' by an external material halo, of the type found in many astrophysical systems, 



o 
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such as elliptical galaxies and globular clusters. Both the Newtonian and the 
general-relativistic dynamics are examined, and in the relativistic case the dy- 
namics of both massive and massless particles are investigated. The halo field 
is given in general by a multipole expansion; we restrict ourselves to multipole 



^ ■ fields of pure order, whose Newtonian potentials are homogeneous polynomi- 

als in cartesian coordinates. A pure n-pole field has n different escapes, one 
of which is chosen by the particle according to its initial conditions. We find 
that the escape has a fractal dependency on the initial conditions for n > 2 
both in the Newtonian and the relativistic cases for massive test particles, but 
with important differences between them. The relativistic motion of massless 
particles, however, was found to be regular for all the fields we could study. 
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The box-counting dimension was used in each case to quantify the sensitivity 
to initial conditions which arises from the fractahty of the escape route. 

I. INTRODUCTION 

It is common in astrophysics to find systems composed by an external material halo 
surrounding an inner region which is either void or has a massive black hole in the center. 
Examples of systems that are realistically described by such a model are elliptical galaxies 
and globular clusters fl]] . The gravitational field in the inner region due to the halo satisfies 
the vacuum field equations, and in Newtonian gravitation it can be expanded in a sum of 
multipole terms of various orders. Due to the intrinsic nonlinearity of general relativity, in 
general this is not true for general relativistic fields. However, in many situations the halo 
can be supposed to be axisymmetric; in this case, the general solution of Einstein's equations 
is known, and the inner relativistic field can also be expanded in a series of multipolar terms. 
In this article, we study the dynamics of test particles in pure multipolar halo fields (without 
any matter in the inner region), both in the Newtonian and relativistic theories, for several 
multipole orders, for material particles and (in the relativistic cases) for light. We find that 
in each case the systems have distinct well-defined escapes. Different escapes are chosen 
by the test particles depending on its initial conditions; if the escape route is a fractal 
function of the initial conditions, the outcome of a particle starting from a given initial 
condition shows a sensitive dependence on these initial conditions. In this sense, we say 
that the system is chaotic @-§|. Associated with this, we have a conveniently defined fractal 
dimension, which quantifies this sensitiveness. Fractal escapes have been found in a number 
of physical systems, including the Henon-Heiles system 0, astrophysical potentials 0, the 
motion of test particles in gravitational waves |M] , classical ionization |3] , motion of particles 



in a fluid ||^, inflationary cosmology [|T0[, and others. We study the fractal dimension for 
different multipole orders and for different energies, in the case of material test particles. 
For simplicity, we confine ourselves to the case of zero angular momentum, when the motion 



is planar. 

This article is organized as follows: in section 2, we investigate the classical multipole 
gravitational field, showing the fractality of the escape basins and calculating the corre- 
sponding dimension for several multipole orders. In section 3, we turn to the relativistic 
case, and investigate the dynamics of both material particles and massless particles; we 
compare the results with those obtained in the classical case. In section 4, we summarize 
our results and draw some conclusions. 

II. NEWTONIAN MULTIPOLE FIELDS 

The Newtonian gravitational field in the inner region due to the mass in the halo satisfies 
the vacuum (Laplace) field equation V^^ = 0, and therefore can be uniquely expanded in a 
sum of multipole terms. If the field is axisymmetric, as we admit from now on, its multipole 
expansion can be written as: 






n=0 

where i?^ = z^ + r^, cos^ = z/R, and r and z are the usual cylindrical coordinates; the Pn 
are Legendre polynomials. If we consider the inner region to be empty of matter, \E' cannot 
have singularities there. This means that bn = for all n in the inner region. Thus: 

oo 

^ = ao + E«"i?"Pn(cos^), (1) 

n=l 

where ao is a constant which can without lack of generality be chosen as zero. 

The multipole expansion of ^ contains only terms that become singular as i? ^ oo (with 
the exception of the constant term ao)- Here we are interested in studying the dynamics of 
the field due to each multipole term separately. We therefore define the n-th order multipole 
field \l'r, as: 



n 



^„ = a„i?"P„(cos^). (2) 



We now consider a test particle moving in the field \Ef„ with no angular momentum in 
the direction of the symmetry axis z; in this case, the particle's motion is restricted to a 
plane that contains the z axis. A convenient set of coordinates in this plane is given by the 
z axis and the extended r axis, which reaches negative values. In order to avoid confusion, 
we shall denote the extended r axis by x, with — oo<x<+oo. \E'is written in terms of x 
simply by replacing it for r. The Hamiltonian of a test particle of unit mass acted upon by 
the field \E'„ is: 

Hn = -(f + z')+mn{x,z). (3) 

The corresponding equations of motion are x = —d'^n/dx, z = —d'^n/dz. Since \&„ is 
a homogeneous potential, satisfying \E'„(Ax, A^) = A"\l'„(x, z), these equations are invariant 
with respect to the similarity transformations x -^ \x, z — * Xz, t — * t/vA, and E -^ X'^E, 
where A is an arbitrary scale factor. Thus, to each orbit in the \E'„(x, z) potential there 
corresponds a similar orbit in the transformed potential ^„(Aa;, Xz), as long as A is positive. 
This means that the dynamics for all positive values of E is the same except for scaling; 
we have to consider only one value of E for each multipole potential \E'„ when we study its 
dynamical properties. We could use this scaling property of the Hamiltonian to make a„ = 1 
in eq. (H), but we prefer to keep the constant arbitrary because of numerical convenience. 

The region in the (x, z)-plane accessible by a particle moving in the potential \E'„ with 
energy E is bounded by the equipotential curves \E'„(x, 2) = E; the particle must satisfy 
\E'„(x, z) < E. This inequality can only be satisfied for a positive energy if \I'„ is positive. 
From equation (H), it follows that the directions 9 satisfying P„(cos^) < are accessible 
for all positive energies, whereas directions satisfying P„(cos^) > 0, if followed far enough, 
reach values of R which violate the inequality \E'n(x, z) < E. Therefore, the directions with 
negative values of P„(cos6') define distinct escapes of the potential \E'„; these escapes are 
separated by inaccessible regions lying at angles with P„(cos^) > 0. From the properties of 
Legendre polynomials, it follows that the potential \E'„ has n distinct escapes. 

As an example. Fig. 1 shows equipotential curves for some energies for the octopole field 



\l/3. The curves for different energies have the same shapes, which is a direct consequence of 
the potential's homogeneity. The three openings in the equipotential curves are seen clearly. 

To each escape we define the corresponding basin, which is the set of initial conditions 
in the phase space that lead to that escape. If the escape chosen by a particle is a fractal 
function of its initial condition, then the basins of the various escapes are mixed in a complex 
fashion, down to arbitrarily small scales, and the final state of the particle (i.e. which escape 
it chooses) can be extremely sensitive to its initial state. This happens when the basin 
boundary is fractal. We then say that the system presents chaos. This is the appropriate 
definition of chaos to open systems, whose motions are not confined to a limited volume of 
the phase space, and this is the definition used throughout this article. 

The directions 6f. which are minima of P„(cos^) correspond to valleys in the potential 
\E'„; the escape regions are centered around these directions. The walls on either side of these 
valleys became increasingly steeper as we go to i? ^ oo. This makes the orbits of escaping 
particles become focused in beams around these directions as they move towards R ^ oo. 
This phenomenon can be seen in fig. 2a, which shows some numerically integrated orbits all 
beginning with initial conditions on x = z = 0, and varying initial velocity direction angle 
0, defined by: 



where v^ and Vz are the initial velocities in the x and z directions, and the velocity modulus 
V is fixed by the conservation of energy: 

V = V2E, 

since ^„(0, 0) = 0. 

The fractal nature of the basin boundary can be observed in a plot of the escape angle 
6e as a function on the initial velocity direction angle. This is shown in Fig. 3a. We had to 
choose a "cutoff distance" Rq such that the integration stops when R reaches Rq. In the case 
of Fig. 3, we chose Rq = 10.0. The discreteness of the values assumed by the escape angle is 



clearly seen in Fig. 3; this is an illustration of the focusing effect. The magnification of a 
tiny region of fig. 3a shown in fig. 3b gives us a convincing proof that the basin boundary 
is fractal. An example of the sensitivity to initial conditions implied by the fractality of 
the basin boundary is given in fig. 2b, where three orbits with differing by less that 
one thousandth of a radian choose different escapes. We see that initially the orbits are 
indistinguishable, and after a time they separate rather abruptly and go to their distinct 
escapes. Figures 2 and 3 were plotted with E = 1 and a^ = 1, but as we explained above, 
the result is exactly the same regardless of the parameters' values. We note that we have 
found no bounded orbits for any of the \E'„ potentials, which might have been expected from 
the shapes of the equipotentials. 

The reason why the escape basins are fractally mixed is that the particle can bounce 
many times off the potential "walls" before escaping. This means that associated with the 
fractality in the basin boundary we have also fractality in the escape time. To see that this 
is so, we plot the time it takes for a particle to escape (without discriminating among the 
various escapes) as a function of its initial velocity angle 0, for a fixed initial position xq, Zq. 
Figure 4 shows the results for the octopole potential, with xq = Zq = 0. The magnification 
shown in fig. 4b shows clearly the fractal nature of the escape time. Notice that fractal and 
non-fractal areas are mixed for all scales. The fractal nature of the escape time is a result of 
the existence of a fractal set of singularities in the escape time function, where the escape 
time diverges. This set of singularities results from the intersection between the set of initial 
conditions (a topological circle in this case) and the set in the phase space composed of 
orbits that never escape for t -^ oo, which in its turn is the stable manifold of the so-called 
invariant set, which is also fractal in our case and is composed of orbits that never escape 
for both t ^ oo and t -^ — oo. 

The sensitivity to initial conditions implied by the fractal nature of the basin boundary 
is made precise with the introduction of the box-counting dimension d, calculated in the 
following way. For a given initial position (which is xq = y^ = in our case), we pick a large 
number of random initial velocity angles 0, with < < 2tt. For each of these we integrate 
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the equations of motion and see to which basin it belongs. We then take two neighboring 
angles (f) — e and + e (with the same initial position), with e being a small number, and 
find to which basin they belong. If all these three initial conditions do not belong to the 
same basin (if they lead to different escapes), then cf) is said to be an "uncertain" angle. The 
number n of uncertain angles found in a sample of N randomly chosen points is a function 
of the displacement e, and for large N the fraction n/N of uncertain points scales as 

I - / « e'-\ (4) 

where d is by definition the box-counting dimension of the intersection of the basin boundary 
with the one-dimensional manifold in which we chose the initial conditions. We numerically 
calculate / for several values of e, and d is directly obtained from the plot of In / versus In e, 
which should be a straight line, and whose inclination a is related to d hy a = 1 — d. 

We have calculated in this way the box-counting dimension for multipole fields with 
orders from 2 to 9. Since the dynamics is the same for every value of the energy, d does not 
depend on E, and is unique for each multipole order. Figure 5 shows the plot of In / versus 
Ine we have obtained for the octopole field (ra = 3). The data fit quite well into a straight 
line, whose inclination (a = 0.77 ± 0.01) gives a dimension d = 0.23 ± 0.01, confirming the 
fractal nature of the basin boundary. Figure 6 shows d as a function of the multipole order n. 
We see that d increases monotonically with n, which means that the fractality of the basin 
boundary increases with n, and thus the dynamics becomes "more chaotic" with multipolar 
fields of higher orders. The dipole and quadrupole fields {n = 1 and n = 2, respectively) 
have regular basin boundaries, and show no chaos in their dynamics. In fact, the potentials 
\l/i and \l/2 are separated using the x,z coordinates, as can be immediately seen from (||), 
which implies that the Hamiltonian is integrable and the basin boundaries are regular. 

A field of the general form (|I]), with the contribution of more than one multipole order, 
is not homogeneous, and its dynamics changes with the energy. We have studied how the 
box-counting dimension changes with E for the field \1' = \E'3 + \E'4, with equal contributions 
of the 3rd and 4th multipole orders. For high values of E, the particle's motion is dominated 



by the 16-pole component of the field, and the equipotentials are practically the same as the 
pure 16-pole field. Accordingly, for these energies the box-counting dimension was found to 
be numerically indistinguishable from the value we obtained for the pure 16-pole field \1'4. 
As the energy is lowered, however, d does not go to its pure octopole field value, as might 
be imagined. In fact, for E = 0.1, we found d = 0.38 ± 0.01; this value of d is higher then 
that of both the pure octopole field and the quadrupole field. Thus, the fractal dimension 
from a field made by the mixture of two multipole orders is not necessarily an intermediate 
value between the dimensions of each pure multipole order separately. 

III. RELATIVISTIC MULTIPOLE FIELDS 

In this section we investigate the dynamics with a general-relativistic generalization of 
the Newtonian multipole fields we studied in the previous section. In general relativity, 
the dynamics results from the space-time metric ds"^ = g^^dx^dx'^ , where g^y is the metric 
tensor and x^ are the space-time coordinates. The equations of motion for a test particle in 
the field g^^ are given by the geodesic equation x^ = F^^x^i;^, where V^j^ is the Christoffel 
symbol, and i;^ = dx'^/dX, with A being a convenient affine parameter. 

For an axially symmetric and static field, the metric is most conveniently written in the 
Weyl form |Tl[] : 



ds' = e^^'dt' - e-2* [e'^ {dr' + dz") + r^dcf?] , (5) 

where r, z and are cylindrical-like coordinates, and t is the time. We use units such that 
c = 1 throughout. The metric functions ^ and 7 depend on r and z only, and satisfy the 
vacuum Einstein field equations: 

Q^2 y, Q^ Q^2 



dj = r 



dr I \dz I 



*--f^- (^) 



Eq. (ID is just Laplace's equation in cylindrical coordinates, and eq. (0) is a quadrature 
whose integrability is automatically guaranteed by eq. (^) . Thus, each axisymmetric solution 
of Newton's field equation corresponds to a general-relativistic metric (^), with 7 being found 
by integration of eq. (|^) with the requirement that the metric be regular in the symmetry 
axis. Based on this, we define the relativistic equivalent to the multipolar fields \l/„ in eq. 
(I) to be the metric (|) with \E' = \E'„ and 7 = 7^ obtained through eq. (0). The integration 
can be performed analytically, and the result is: 

fl r -1 

7„ = -a-'R^- [PnicosO) - PLiicosO)] , (8) 

with i?^ = r^ + z^. Since there are many possible relativistic space-times which in the weak 
field limit correspond to a given Newtonian field, the choice of the metric is not unique. We 
have chosen the metric defined by \E'„ and 7„ because it is the simplest choice, and because 
the field obtained in this way shares important features with the Newtonian field, as we 
shall see shortly. 

As in the Newtonian case, the invariance of the system with respect to time translations 
and axial rotations implies the existence of two constants of motion, which can be identified 
as the energy E and the ^;-component of the angular momentum L^. They can be found by 
using the Lagrangian associated with the metric (1): 



2L = e^^P - e-^^ 



e^^ 



(r^ + i2) + ^2^2j ^ ^g^ 



where the dot denotes differentiation with respect to the affine parameter. Since both \E' and 
7 depend only on r and z, the Lagrangian L does not depend on t and 0; they are "cyclical 
coordinates" in the language of Hamiltonian dynamics, and they are associated with the 
conserved quantities: 

^ - t - e», (10) 

i. = g = -r'e-«4.. (11) 



Since the Lagrangian function (j^) has only kinetic terms, and since L does not depend 
exphcitly on the affine parameter, the Hamiltonian of this system is simply L itself, and its 
conservation gives us yet another constant of motion: 

where niQ is the test particle's rest mass. In the case of light, we have rriQ = 0; for massive 
particles, we can without loss of generality make mo = 1. Rewriting L with the help of eqs. 
(0) and (inD, we have: 

E2e-2* - L^r-^e^* - e2(^-*) (f^ + ^2^ = ^^ (12) 

with 6 = for particles with zero rest mass, and 6=1 for massive particles. As in 
the Newtonian case in the previous section, we will restrict ourselves to motion with zero 
angular momentum (L^ = 0), when the motion is planar. In complete analogy to what we 
did for the Newtonian potential, we define a new coordinate x, which is no more than r 
extended to negative values; the motion of our zero angular momentum test particle occurs 
on the (x, z)-plane. All the expressions above are written in terms of x simply by replacing 
r by X everywhere. 

The accessible region of the (x, 2;) -plane is found by setting a; = i = in (0); taking 
into account that L^ = 0, we have: 

E\x = z = 0) = U{x, z) = 6e^^. (13) 

The accessible region for a particle with energy E is given by U{x, z) < E"^. For massless 
particles {6 = 0), U = 0, and this condition is always satisfied; thus, the whole (x, 2;)-plane 
is accessible to null geodesies. This is not so for massive particles {6 = 1); the regions they 
are able to reach on the xz plane depend on their energy and are determined by \&. 

We now consider the pure multipole potentials \&n = ctnR"' Pn{cos 9) , with 7„ given by 
(^). We first treat the case of massive particles; the motion of light will be studied later. The 
effective potential is U = e^* = gu, where gu is the time-time component of the metric. It is 
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clear that the dynamics in the relativistic field {"^n, In) is more complex than its Newtonian 
counterpart. In particular, since neither the effective potential (p!3| ) nor 7„ are homogeneous 
functions, and furthermore the relativistic equations of motion depend on the velocities x'^, 
the orbits in the relativistic field do not satisfy the homogeneity property that guarantees for 
the orbits in the Newtonian field that the dynamics is the same for all energies, up to a global 
scaling. We thus expect the dynamics to change with the energy in the relativistic case; this 
is a fundamental difference between the two theories. In spite of this, we see from eq. ( |13D 
that the curves of constant effective potential are still given by \E'„ = const (although the 
constant is no longer just the energy). This means that the escapes in the potentials \&„ 
are the same as those in the Newtonian fields. As -R ^ oo along the asymptotically allowed 
directions, given by P„(cos^) < (exactly as in the Newtonian case), we have U = gu ^ 0, 
that is, we have an infinite red shift (as in section 2, we assume a„ > 0): particles are 
attracted to these regions. On the other hand, for asymptotically forbidden regions, with 
Pn{cos9) < 0, we have gu -^ oo, an infinite blue shift: particles are attracted away from 
these regions. 

We now study the dynamics in this field using the same techniques we used in the 
Newtonian case. We integrate numerically orbits with initial position xq = zq = and 
variable velocity angle 0, with the magnitude v of the velocity being fixed by the conservation 
of energy. Since \E'n(0, 0) = 7„(0, 0) = for all n, from eq. (|12D we see that: 



Xq = V sm (p Zq = V cos ( 



with 



t,2 = ^2 _ ^^ 



The equations of motion for the independent variables x and z are obtained directly from 
the geodesci equations: 

1 



^^-^ 



'g%E' + U {i' - i') + V,zii] ; (14) 
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z = -^ [g%E^ + /,. (i' - x^) + 2/,,xi] , (15) 

with / = e^*^*"''') and (7** = e~^*; these equations are vahd for L^ = 0. Together with the 
quadrature (0), the above equations give the complete dynamics of the system. 

Similar to what we found in the previous section, the orbits are focused along directions 
corresponding to the minima of P„(cos 6). Plots of the escape angle for n > 2 are very similar 
to Fig. 3, showing clearly that the escape has a fractal dependency on the initial conditions. 
As in the Newtonian case, the dipole and quadrupole fields (n = 1 and 2, respectively) 
show a regular basin boundary. However, the equations of motion are not separated in x, z 
coordinates, as is the case for the Newtonian fields. 

As we explained above, we expect the fractal basin boundary present in the n > 2 fields 
to change with the energy in the relativistic case, as opposed to the Newtonian one. As in 
the Newtonian case, we have found no bounded orbits, neither for massive particles nor for 
massless ones. 

We calculate the fractal dimension of this set in the same way we did previously for 
the motion in the Newtonian field. In that case the dimension did not depend on the 
energy, because the dynamics is the same for all (positive) energies. This is not true for 
the relativistic dynamics, and we expect the dimension to change with the energy. This is 
indeed what we have found. As an example, we take the field (\E'4, 74). We consider particles 
with energies higher than 1, corresponding to positive energies in the Newtonian case. For 
E close to 1 (corresponding to energies close to in the Newtonian field), we find that the 
box-counting dimension is equal to its Newtonian value. For instance, for E = 1.01, we find 
d = 0.33 ±0.02, which is (within the error) the same value we found for the Newtonian field 
\E'4. For higher energies, however, d changes with E: for E = 1.5, we have d = 0.30 ± 0.02, 
and for E = 2.0, d = 0.26 ± 0.02. We have calculated d for more values of E than we 
show here, and it seems that d decreases monotonically with E, and the system becomes 
more chaotic for large energies; this is the contrary to what is usually the case in bounded 
dynamics, where commonly the chaos grows with the energy. The same behavior is observed 
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in the octopole field (n = 3), but in this case the variations in the dimension are smaller. 
Because of the complexity of the equations of motion for high multipole orders, we have 
not been able to calculate d with enough accuracy for n larger than 4, but we think this 
behavior of d might be a general feature of the dynamics for the \l/„ fields. 

In the case of null geodesies, the equations of motion are again given by eqs. (|Tip and 
(p!5|), but the initial conditions are required to satisfy the constraint ([T^ ) with 5 = (and 
L2 = in our case). This allows us to scale the affine parameter so as to eliminate E from 
the equations. This means that the planar motion of light depends only on the chosen field 
\l/„, and not on any other parameter, unlike the motion of massive particles. 

Even though the whole (x, 2;)-plane is accessible to massless particles, as we have seen, 
null geodesies also show focusing effects similar to that shown by massive particles. An 
example is given in Fig. 7, where the escape angle is plotted as a function of the initial 
velocity angle for null geodesies, for the octopole field n = 3. We see clearly that the values 
assumed by the escape angle are indeed discrete. We do not see, however, any region where 
the escape angle shows rapid variations, that is typical of fractal functions. Magnifications 
of several regions of Fig. 7 show no small-scale structures, indicating that the escape angle 
is not fractal. This is confirmed by direct calculation of the box-counting dimension, which 
gives d = 1 io within statistical accuracy. This result also holds for n = 4 and n = 5, which 
is the highest multipole order for which we have been able to obtain reliable results. We 
speculate that this result might be true for all n. 

IV. CONCLUSIONS 

We have studied the dynamics of test particles in pure multipole fields generated by 
external halos of matter, when the particles have zero angular momentum and the motion 
happens on a plane. Such a system could be used to approximate observed astronomical 
fields, such as the interior of galaxies. We have considered both the Newtonian and the 
general-relativistic fields, and we have compared the results of the two theories. 
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For the Newtonian fields, the homogeneity of the pure potentials "^n implies that the 
dynamics is independent on the energy. We have found that the escape chosen by the test 
particles is a fractal function of the initial conditions for all pure multipole fields of order 
higher than 3. For n = 1,2, the Hamiltonian is separable and the basin boundaries are 
regular. We have calculated numerically the box-counting dimension d associated with the 
escape function. Since the dynamics is the same for all energies, d depends only on the 
multipole order n. We have found that d is an increasing function of n for n up to 9; higher 
values of n could not be studied because of numerical limitations. The increasing in the 
intensity of chaos with n can be understood by the fact that fields of higher order have 
a higher number of escapes. We have also studied the dynamics of a field made by the 
superposition of two different multipole orders, specifically, the field \1' = ^3 + \E'4. For this 
field, the dynamics depends on the energy, and we found that, although for high energies 
the box-counting dimension tends to the value of the \E'4 field, for low energies d assumes 
values higher than the dimensions of both ^3 and \E'4, maybe somewhat surprisingly. 

We now discuss the results found for the relativistic fields, starting by the dynamics of 
massive particles. The relativistic pure multipole fields have not the homogeneous property 
shown by their Newtonian counterparts, and their dynamics does depend on the energy, and 
so does the escape function. We found that the escape chosen by the particle is a fractal 
function of the initial conditions for fields of order higher than 2 also in the relativistic case. 
The dipole and quadrupole fields have regular basin boundaries, but unlike the Newtonian 
fields, the relativistic equations of motion are not directly separable. For n > 2, we expect 
the box-counting dimension d to change with the energy. We have calculated d as a function 
of the energy, and we have found that for low energies, d goes to the Newtonian values, as 
might have been expected. As E increases, d was found to decrease monotonically, for all 
the energies for which we have been able to obtain reliable results. We have obtained these 
results for n = 3, 4 and 5, and this seems to be a general feature of the dynamics of these 
fields. 

For massless test particles, the escape function was found to be regular (non-fractal) for 
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n up to 5; at least as far as the fifth multipole order, null geodesies are not chaotic. 
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FIGURE CAPTIONS 

Fig. 1 Equipotential curves for the Newtonian octopole field \l/3. From tlie inside out, the 
energies are E = 0.25, E = 1.0 and E = 4.0. 

Fig. 2 Orbits in the Newtonian field \1'3 all starting from x = z = 0, with different initial 
directions and E = 1.0, a^ = 1. Fig. 2b shows three orbits with initial velocity angles 
differing by 1/5000 radian; the circle marks the starting point x = z = 0. 

Fig. 3 Escape angle 6e as a function of the initial velocity angle for the octopole field. 
Fig. 3b is a magnification of a small region of Fig. 3a. 

Fig. 4 Escape time as a function of the initial velocity angle for the octopole field. Fig. 
4b is a magnification of a small region of Fig. 4a. 

Fig. 5 In / vs. In e for the Newtonian octopole field. The straight line is the result of the 
fitting. 

Fig. 6 Box-counting dimension d for Newtonian multipole field of order n. The error bars 
come from the statistical errors in the calculation of d. 

Fig. 7 Escape angle 6e as a function of the initial velocity angle for the null geodesies in 
the relativistic octopole field. 
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